Is Multiscaling an Artifact in the Stochastically Forced Burgers Equation? 
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We study turbulence in the one-dimensional Burgers equation with a white-in-time, Gaussian 
random force that has a Fourier-space spectrum ^ 1/fc, where k is the wave number. From very-high- 
resolution numerical simulations, in the limit of vanishing viscosity, we find evidence for multiscaling 
of velocity structure functions which cannot be falsified by standard tests. We find a new artifact in 
which logarithmic corrections can appear disguised as anomalous scaling and conclude that bifractal 
scaling is likely. 
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Homogeneous, isotropic fluid turbulence is often char- 
acterized by the order-p velocity structure functions 
5p(£) = {[{v{x + i) - v{x)} ■ (f)]P), where v{x) is the 
velocity at the point x and the angular brackets de- 
note an average over the statistical steady state of the 
turbulent fluid. For separations (, in the inertial range, 
rjd -^ £ -^ L, one has Sp{£) ~ 1''^. Here rjd is the small 
length scale at vi^hich dissipation becomes important; L is 
the large length scale at which energy is fed into the fluid. 
The 1941 theory (K41) of Kolmogorov [l| predicts simple 
scaling with exponents Cp'^^ = p/3- By contrast, exper- 
iments and direct numerical simulations (DNS) suggest 
multiscaling with ^p a nonlinear, monotonically increas- 
ing, convex function of p, not predictable by dimensional 
analysis ^]. However, the Reynolds numbers achieved in 
DNS are limited, so the exponents (p have to be extracted 
from numerical fits over inertial ranges that extend, at 
best, over a decade in i. The processing of experimental 
data - although they can achieve much higher Reynolds 
numbers - involves other well-known difficulties [3| . It is 
important therefore to establish, or rule out, multiscal- 
ing of structure functions in simpler forms of turbulence, 
such as passive-scalar, passive-vector or Burgers turbu- 
lence. Significant progress, both analytical and numeri- 
cal, has been made in confirming multiscaling in passive- 
scalar and passive- vector problems (see, e.g., Ref. 4] for 
a review) . The linearity of the passive-scalar and passive- 
vector equations is a crucial ingredient of these studies, 
so it is not clear how they can be generalized to fiuid 
turbulence and the Navier-Stokes equation. 



Here we revisit the one-dimensional. Burgers equa- 
tion with stochastic self-similar forcing, studied earlier 
in Refs. |a,l3- It is by far the simplest nonlinear partial 
differential equation (PDE) that has the potential to dis- 
play multiscaling of velocity structure functions [a| ; and 
it is akin to the Navier-Stokes equation. In particular, 
we investigate the statistical properties of the solutions 



to 



dtu + udxU = vdxxU + f{x,t), 



(1) 



in the limit of vanishing viscosity j/ — > 0. Here u is the 
velocity, and f{x,t) is a zero- mean, space-periodic Gaus- 
sian random force with 

(/>i, ti)/>2, ^2)) = 2Do\k\^d{h - t2)S{h + h) (2) 

and f{k, t) the spatial Fourier transform of /(x, t). We re- 
strict ourselves to the case /? = — 1 and assume spatial pe- 
riodicity of period L. Earlier studies |a,la| suggested that 
Eqs. 1^ and (0), with /3 = — 1, show a nonequilibrium 
statistical steady state with bifractal scaling: this means 
that velocity structure functions of order p < 3 exhibit 
self-similar scaling with exponents p/3 and implies a K41- 
type —5/3 energy spectrum, predictable by dimensional 
analysis, whereas those of order p > 3 have exponents all 
equal to unity being dominated by the finite number of 
shocks, with 0{L^^"^) strength, typically present in the 
periodic domain; this bifractal scaling is somewhat simi- 
lar to that observed when the Burgers equation is forced 
only at large spatial scales |3,|a|- 

We overcome the limitations of these earlier stud- 
ies la 01 by adapting the algorithm of Refs. [3, ll^ to 
develop a state-of-the-art technique for the numerical so- 
lution of Eqs. (Q) and 10), in the p —^ limit. This yields 
velocity profiles (Fig.^a) with shocks at all length scales 
resolved. Structure functions [Figs. I^b) and iQc)] ex- 
hibit power-law behavior over nearly three decades of r; 
this is more than two decades better than in Ref. |5^ . In 
principle it should then be possible to measure the scaling 
exponents [Figs. (Pb)] with enough accuracy to decide 
between bifractality and multiscaling. A naive analysis 
[Fig. |(21a)] does suggest multiscaling .I^]. However, given 
that simple scaling or bifractal scaling can sometimes be 
mistaken for multiscaling in a variety of models [l^llSjj it 
behooves us to check if this is the case here. We describe 
below our numerical procedure and the various tests we 
have carried out. 
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FIG. 1: (a) Representative snapshots of the force / and the velocity u (jagged Une), in the statistically stationary regime; the 
velocity develops small-scale fluctuations much stronger than those present in the force. Log-log plots of the structure function 
Sp'"^ (r) versus r for A'' = 2^" and (b) p = 3 and (c) p = 4. The straight line indicates the least-squares fit to the range of scales 
limited by the two vertical dashed lines in the plots. The resulting multiscaling exponents ^p (see text) are shown by horizontal 
lines in the insets with plots of the local slopes versus r. 



In our simulations we use L = 27r and Dq = 1 without 
loss of generality. The spatial mesh size is 8x = L/7V 
where 'N is the number of grid points. To approximate 
the forcing, we use the "kicking" strategy of Ref. |9j 
in which the white-in-time force is approached by shot 
noise. Between successive kicks we evolve the velocity 
by using the following well-known result on the solutions 
to the unforced Burgers equation in the limit of vanish- 
ing viscosity (see, e.g., Ref. Q): the velocity potential '0 
(such that u = ^dxip) obeys the maximum principle 



4^{x, t') — max 
y 



^iy,t)- 



jx - yf 
2 {f - 1) 



t' > t. 



(3) 



The search for the maxima in Eq. ^ requires only 
0(iVlog2 A^) operations [lOj because, under Burgers dy- 
namics, colliding Lagrangian particles form shocks and 
do not cross each other. At small scales we want to have 
unforced Burgers dynamics with well-identifiable shocks. 
At least four mesh points are needed for unambiguous 
identification of a shock; since the maximum wavenum- 
ber is N/2, we set f{k,t) = beyond an ultra-violet 
cutoff A = N/8. 

Specifically, at time i„ — nSt we add fn{x)\/6t to 
the Burgers velocity u{x,t), where the /„(a;)s are in- 
dependent Gaussian random functions with zero mean 
and a Fourier-space spectrum ~ 1/fc, for fc < A. The 
time step 5t is chosen to satisfy the following conditions 
{6x/2uo) < St ~ {l/LA)^/^{L/uo), where Uo is the char- 
acteristic velocity difference at large length scales 0{L). 
The first ensures that a typical Lagrangian particle moves 
at least half the mesh-size in time St (otherwise, it would 
stay put)jl^]; the second, which expresses that St is com- 
parable to the turnover times at the scale L/A, guaran- 
tees that, at scales larger than L/A the time-step, St, 
is small compared to all dynamically significant times, 
but still permits the formation of individual shocks at 
smaller scales. Finally, as our simulations are very long, 
for the stochastic force we use a good-quality random- 



Run 


N 


St 


A 


TL 


Ttr 


J- av 


Bl 
B2 
B3 


020 
2l8 
216 


5 X 10-* 
1 X 10-* 
1 X IQ-* 


2'' 

2l5 
ol3 


1.0 
1.0 
1.0 


2.0 
2.0 
2.0 


22 

20 
120 



TABLE L Different parameters used in our runs Bl, B2 and 
B3. TL = L/uo is the equivalent of the large-eddy-turnover 
time. Data from Ttr time steps are discarded so that tran- 
sients die down. We then average our data over a time Tav 



number generator with a long repeat period of 2™ due 
to Knuth |ll| . The main characterisitics of the runs per- 
formed are summarized in Table Q] 

In addition to the usual structure functions, we have 
used S^^^ir), defined by 



Su{x,r) = u{x + r) — u{x). 



ip 



(4) 
(5) 



from which we extract the exponents S,p- For each value 
of N we have calculated £,p for p — ra/A, with integers 
\ < m < 20. Figure (|21 a) summarizes the results of 
our calculations concerning the exponents ^p, for N = 
2^^, 2"'^*, and 2^°; any systematic change in the values of 
these exponents with N is much less than the error bars 
determined by the procedure described below. Thus in 
all other plots we present data from our simulations with 
N — 2^^ grid points. The representative log-log plots of 
Figs. Qb) and (Pc) of S^'""(r) for p = 3 and 4 show 
power-law regimes that extend over nearly three decades 
of r/L. We obtain our estimates for the exponents ^p as 
follows: for a given value of p we first determine the local 
slopes of the plot of log S^^ versus log r by least-squares 
fits to all triplets of consecutive points deep inside the 
power-law regime |l7| . These regions extend over one 
and half decade of r/L as shown in Figs. iQb) and l^c). 
The value of ^p we quote [Fig. |(21a)] is the mean of these 



local slopes; and the error bars shown are the maximum 
and minimum local slopes in these regions. 

Figure (|21 a) shows that our results for ^p, indicated 
by circles for N = 2"^°, deviate significantly from the 
bifractal-scaling prediction (full lines). As we shall see 
this deviation need not necessarily imply "multiscaling" 
for the structure functions. We have also considered the 
possibility of already well-understood artifacts, such as 
the role of temporal transients [l3| and finite-size effects 
which can round sharp bifractal transitions [l3| , and de- 
cided that they do not play any major role in the present 
problem |ia |. 

Consider P'^{s), the cumulative probability distribu- 
tion function of shock strengths s. Simple scaling argu- 
ments predict P'^'{s) ^ s'* , with 7 = —3, which follows 
by demanding that P'^(s) remain invariant if lengths are 
scaled by a factor A and velocities by A~^/'^. One of 
the signatures of multiscaling would be deviations of 7 
from this scaling value. However, the following argu- 
ment favors 7 = — 3 : the total input energy, E\n = 
Ji^ D{k)dk ^ In A, where A is the ultra-violet cutoff. 
In the limit of vanishing viscosity, the energy dissipa- 
tion in the Burgers equation occurs only at the shocks 
and is proportional to the cube of the shock strength [g , 
so the total energy dissipation is fJ ^ J /"'^'' P{s)s'^ds. 
Here P{s) = {dP'^{s)/ds) ~ s'^~'^ is the probability dis- 
tribution function (PDF) of shock strengths and Smin and 
Smax are, respectively, the minimum and maximum shock 
strengths. A steady state can occur only if Ein and fi 
have the same asymptotic properties as A ^ 00, i.e., 
H. ^ ln(A); this requires 7 — —3. By contrast we find 
7 ~ —2.7 [Fig. Elb)] by a naive least-squares fit to the 
tail of P'^{s) 19] . This suggests that the results of our 
simulation are far from the limit A — > 00, although more 
than a million grid points are used. Hence the "anoma- 
lous" exponents in Fig. |(5]a) might well be suspect. 

To explore this further, consider the third-order struc- 
ture function of velocity differences, without the absolute 
value, namely, 5*3 (r) = (Su^). From Eqs. iQJ and (jSJ fol- 
lows the exact relation 



Is^ir) = J^ Fiy)dy, 



(6) 



where F{y) is the spatial part of the force correlation 
function, defined by {f{x + y,t')f[x,t)) = F{y)5{t - t'). 
We obtain this analog of the von Karman-Howarth rela- 
tion in fluid turbulence by a simple generalization of the 
proof given in Ref. 3>] for the Burgers equation forced de- 
tcrministically at large spatial scales. An explicit check 
of Eq. @ provides a stringent test of our simulations 
[Fig. (0 c) inset]. Furthermore, Eq. 10) implies that 
S^{r) ^ rlog(r) for small r and thus should display 
significant curvature in a log-log plot, as is indeed seen 
in Fig. (0 c). By contrast Sf'^r) [Fig.(j3 c)] displays 
much less curvature and can be fitted over nearly three 
decades in [r / L) to a power law with an "anomalous" 



exponent of 0.85. This anomalous behavior is prob- 
ably an artifact as we now show. Let us define the 
positive (resp., negative) part of the velocity increment 
5^u (resp., 5~u) equal to 5u when ^u > and to zero 
when ^u < (resp., to 5u when 5u < Q and to zero 
when 5u > 0). Obviously Sair) = {{S+u)'^) + {{S-uf), 
whereas Sf%r) = {{d+uf) - (((5"u)3). The log-log plot 
of (((5+u)3) = {l/2)[S3{r) + Sl'°''{r)] in Fig. ^c) is much 
straighter than those for 6*3 (r) and Sf^^{r) and leads to 
a scaling exponent 1.07 ± 0.02, very close to unity. For 
a moment assume that {{d~^u)^) indeed has a scaling ex- 
ponent of unity. Given that Ssir) has, undoubtedly, a 
logarithmic correction, it follows that S^^^{r) has (except 
for a change in sign) the same logarithmic correction in 
its leading term (for small r) but differs by a subleading 
correction proportional to r. This subleading correction 
is equivalent to replacing r log(r) by r log(Ar) for a suit- 
ably chosen factor A. In a log-log plot this shifts the 
graph away from where it is most curved and thus makes 
it straighter, albeit with a (local) slope which is not unity. 

An independent check of {{6~^u)^) ^ r is obtained 
by plotting the cumulative probabilities, 'I''^, of posi- 
tive and negative velocity increments (for a separation 
r = 800(5x) in Fig. ^h) |20j. For positive increments 
^"^ falls off faster than any negative power of Su, but, 
for negative ones, there is a range of increments over 
which ^"^ ^ \Su\~^, the same —3 law seen in P'^(s) ear- 
lier. Indeed the negative increments are dominated by 
the contribution from shocks. Just as P'^{s) has cut- 
offs Sinin and Sinax, ^'^ has cutoffs u^-^{r) and Wmax('^) 
for negative velocity increments. Since $'^ falls off as 
[(5u]~^, u~ax <3an be taken to be c»; furthermore as the 
PDF of velocity differences, $(^u) = d{^'^) / d{5u) , must 
be normalizable, we find u^^^{r) ^ r^'^. We now know 
enough about the form of $ to obtain, in agreement with 
our arguments above, that 6*3 (r) w — Ar ln(r) -I- Br and 
S'f'='^(r) w Arln(r) + Br, whence {(S+u)^) « Br. The 
presence of this cutoff yields a logarithmic term in both 
53 and S^^^ but with different sign agreeing with the 
arguments given in the previous paragraph. 

By a similar approach, we find S4{r) « Cr — Dr"^'^, 
where C and D are two positive constants. The negative 
sign before the sub- leading term (r^'^) is crucial. It im- 
plies that, for any finite r, a naive power-law fit to S4 can 
yield a scaling exponent less than unity. The presence of 
sub-leading, power-law terms with opposite signs also ex- 
plains the small apparent "anomalous" scaling behavior 
observed for other values of p in our simulations. A simi- 
lar artifact involving two competing power-laws has been 
described in Ref. [ij]. 

In conclusion, we have performed very-high-resolution 
numerical simulations of the stochastically forced Burg- 
ers equation with a 1/k forcing spectrum. A naive inter- 
pretation of our data shows apparent multiscaling phe- 
nomenon. But our detailed analysis has identified a 
hitherto- unknown numerical artifact by which simple bis- 
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FIG. 2: (a) The multiscaling exponents ^p versus order p for Eqs. Q and J^J with A'' = 2^^(o), 2^®(*), and 2^''(o) grid points. 
Error bars (see text) are shown for the case A'^ = 2^". The deviation of ^p from the exponents for bifractal scahng (full lines), 
shown as an inset, suggest naive multiscaling. (b) Log-log plots of the cumulative probability distribution function P'^{s) versus 
shock strengths s obtained from an average over 1000 snapshots. A least-squares fit to the form P''{s) ~ s"', for the dark points 
in the range —5 < log j^g [-PC'S)! S: —2.5, yields 7 — —2.70; the simple-scaling prediction 7 — —3 is indicated by the straight 
line. The inset shows log-log plots of the cumulative probability distribution function, $'^[(5it(r)], (dashed line : positive Su, 
continuous line : negative 5u) versus the velocity difference Su{r) for length scale r — 800Sx. (c) Log- log plots of 5*3 (r) (crosses), 
S3°^{r) (dashed line) and {{S'^u)^) (squares) versus r. The continuous line is a least-square fit to the range of points limited by 
two vertical dashed lines in the plot. Inset: An explicit check of Eq. JSJ from our simulations, plotted on a log-log scale. The 
dashed line is the right-hand side of Eq. Q; the left-hand side of this equation has been obtained for A'^ — 2^" (o) (run Bl). 



caling can masquerade as multiscaling. Our work illus- 
trates that the elucidation of multiscaUng in spatially 
extended nonlinear systems, including the Navier-Stokes 
equation, requires considerable theoretical insight that 
must supplement state-of-the-art numerical simulations 
and experiments. 
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